Simulations of Pattern Formation in Vibrated Granular Media 
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We present simulations of peak pattern formation in vibrated two- 
dimensional (2D) granulates and measure the dispersion relation of the pat- 
\ tern for various frequencies, accelerations, cell sizes, and layer heights. We 

report the first quantitative data from numerical simulations showing an in- 
teresting dependence of the pattern wavelength on the acceleration and the 
system size. Our results are related to recent experimental findings and the- 



oretical predictions for gravity waves. 
Pacs: 46.10. +z, 47.20-k, 05.60+w 
The dynamical properties of non-cohesive granular media have attracted a lot of interest 
in recent years [|TJ. Vibrated granular assemblies show a variety of possible responses like 
surface fluidization [0-0], convection J7J, heaping |J, and surface waves P~P1) an phenomena 
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being observed in both 2D and 3D systems. In 3D experiments on vibrated layers of sand, 
Melo et al. [11]] find surface patterns similar to the surface waves obtained by parametric 



excitation in regular fluids, i.e. the Faraday instability ||15|| . For a review concerning gravity 
waves in fluids see Ref. |lf| and refs. therein. The patterns in granular materials, viewed from 



the top of the 3D cell, display regular structures such as stripes, squares or hexagons. Recent 



experiments in a reduced 2D geometry on layers of aluminum beads [12] have also shown 



the formation of a peak pattern instability. From these measurements of the parametric 
excitation, a dispersion relation was reported, analoguous to the findings of Melo et al. 



TT| . No clear mechanism for the instability has been given so far, but results seem to 
indicate a behavior specific to granular assemblies. At high frequencies or great layer heights, 
the wavelength saturates at a value independent of the excitation frequency. In a more 



recent study a scaling was proposed for the high frequency limit |T7|]; the argument accounts 
for a dissipation mechanism due to some granular viscosity but no clear evidence for the 
validity of the scaling was given either. Here, we report simulations of vibrated 2D arrays of 



polydisperse spheres and compare with experiments ||12|| . Using material parameters close 
to experimentally reported values, we observe that the instabilities appear rapidly and that 
the dispersion relation for different geometries, excitation, material parameters, and initial 
conditions can be studied in detail. 

The model system consists of iV spheres of diameter d{ (i = 1,...,N) randomly chosen 
from the interval [d(l — w), d(l + w)] with d = 1.5 mm and w = 0.1. Testing different values 
of w we observed the surface waves even in the extreme monodisperse case w — 0. We use 
cells of horizontal width L and not limited in height. The container moves with a vertical 
trajectory: z(t) = Asin(ujt), where A is the amplitude and / = u/(2tt) the frequency of 
excitation. The maximum acceleration of the bottom plate is defined as the dimensionless 
quantity T = Au 2 /g, with the gravitational acceleration g. The layer thickness is H = Nd/L 
with the dimensionless system width L/d. 
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Here, an event driven (ED) method is used to simulate the dynamics of rigid h ard spheres 
with no intrinsic material elasticity, i.e. the duration of a contact is zero (t c = 0). This choice 
is made for two reasons: Firstly, it is crucial to show that the instability presented here, is 
decoupled from a possible parametric excitation of collective elastic modes which might be 
generated in the elastic network formed by a dense packing of soft spheres (t c > 0), i.e. the 



"detachment effect" [0. Secondly, in the range of response where the typical separation 



between the beads is not too small, the ED simulation scheme is quite efficient. The method's 
principles are the following: under the influence of gravity, the particles follow a parabolic 
trajectory, until an event occurs. An event is either a collision between two particles, or 
a collision of one particle with a wall or the bottom plate. From the velocities just before 
the event, the velocities after this event are computed, accounting for the energy loss due 
to friction, and some inelasticity of the material. In the tangential direction, we account 
for friction, using the coefficient of friction fx and the maximum tangential restitution (3 . 
This simplified description of tangential dissipation introduces a coupling between the linear 
and the rotational degrees of freedom || and is consistent with recent experimental results 
on colliding particles JT9|. Simulations performed with different fx and (3q values show the 
instability, even for no rotational coupling at all, i.e. fx = 0. This means that the rotational 
degree of freedom of the particles is not crucial for the instability to occur. Furthermore, the 
patterns also occur when dissipation and friction at the walls is switched off, proving that 
the pattern forming instability is not influenced by the wall's properties, as e.g. convection 
is. Switching off friction with the bottom leads to less stable patterns in the sense, that the 
peaks move in the horizontal direction more easily. 



In contrast to previous ED simulations of granular assemblies P,[5|, |20|j21|] , we implement 
a dissipation model which uses a velocity dependent restitution coefficient. Such a model 
is qualitatively consistent with experimental measurements of binary collisions, reporting 



a restitution coefficient approaching unity with decreasing velocity |f22|-[24|j. The specific 
model we consider here is the limiting case of a visco-elastic interaction law for the con- 
tact of spherical particles. The variation of the surface of contact during the interparticle 



penetration causes a non-linear elastic force, i.e. the Hertz model [25], and a non-linear 



dissipative force, i.e. the Kuwabara-Kono model |23[. Thus we use here the coefficient of 
normal restitution e(u) = 1 — EqIu/uo) 1 , with the relative velocity in the normal direction u, 
and the power 7 = +1/5. In order to model aluminum spheres we use Eo = 0.4, which leads 
for typical velocities of uq = lm/s to e — 0.6. Larger velocities yield smaller coefficients 
of restitution. Other models, e.g. based on the hypothesis of plastic deformation, lead to 
qualitatively similar behavior f26f , but so far, experiments do not provide a discriminant 
conclusion about the exact velocity dependence |T^J2^J2^|. 

However, this systematic decrease of dissipation with decreasing average relative energy 
is not sufficient to keep the collision frequency small under all conditions. The divergence of 
collision frequency, connected to a (possibly local) loss of relative energy is usually referred 
to as the "inelastic collapse". In order to keep our simulations out of the regime of the 
inelastic collapse, we also introduce a cut-off time t x , prohibiting dissipation for a second 
collision within a time-interval shorter than t x . In order to test the sensitivity of the method 
to the exact value of this cut-off time, we have compared identical simulations using only 
different t x values t x = 0s, 10~ 6 s, 10 _5 s, and 10 _4 s. Except for the largest t x values, we 
get quantitatively the same results. Note that surface waves may occur also in the case of 
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traditional ED simulations with extreme values t x = Os and e = const, i.e. 7 = 0. In general, 
the computational effort decreases with increasing cut-off time, and allows simulations with 
rather large values of H. The authors are aware that more detailed studies are neccessary 
to investigate all the implications of the ED-extensions used here; however, this problem is 
out of the scope of this paper. 

In Fig. [I] snapshots of a typical simulation with N = 600 particles, in a box of width L/d 
= 100, vibrated with / = 10 Hz and the acceleration T = 3.6 are plotted. The parameters 
of dissipation and friction are here e = 0.4, e 0w = 0.2, fj, — [i w — 0.2, and f3 = [3 0w = 
0.0, where the index w indicates the particle-wall interaction parameters, corresponding to 
strong particle-particle and weaker particle-wall dissipation. In Fig. [I], we present a time- 
series ranging from t = 1.30s to 1.48s. We observe, like in the experiment, a parametric 
response of the layer with a period 2T = 2/ /. When the bottom moves up, see (a) and (e), 
the array is compressed and the peaks vanish, see (b) and (f). The array separates from the 
bottom plate after the latter accelerates downwards, see (c) and (g), and the peaks grow, 
see (d) and (h), until the array hits the bottom again. Note that the position of the peaks 
is interchanged with the position of the valleys from one period to the next. An arch-like 
structure below the array, just before the collision with the bottom plate, is also visible, see 
(h). This behavior is in agreement with the experimental findings of Clement et al. JT2{]. 



(a.)t- 1.30 s 



(b)t= 1.32 s 



(c)t= 1.34 s 



(d)t= 1.38 s 



(e)t= 1.40 s 



(f)t= 1.42 s 



(g)t= 1.44 s 



ISP*? 



(h)t = -1.48 s* 



FIG. 1. Snapshots of a typical simulation over two periods from t = 1.30 s to 1.48 s. The 
parameters are N = 600, L/d = 100, / = 10 Hz, T = 3.6, Eq = 0.4, sow = 0.2, \i = 0.2, and (3q = 
0. The dashed line indicates z = 0. 



To perform systematic quantitative measurements on the wavelength of the observed 
pattern, we monitor the behavior of the horizontal particle-particle correlation function: 
C x>x {x) = ( L _^)jv2 ^iLi^jLi 8 [x — \xi — Xj\], with the delta function S[x] = 1 for x = and 
8[x] =0 elsewhere. If the instability is present, this function displays some modulation in 
x with a first maximum which we indentify with the wavelength L x (t) of the surface waves. 
The modulation is strongest just before the collision of the granular layer with the bottom 
plate. We trace L x (t) over 20 to 50 periods and get L x , the averaged wavelength. 

Clement et al. |12| observe from experiment little influence of the acceleration T on 
the dispersion relation of the waves. This was found for aluminum beads in a range of 
accelerations where the frequency of layer-bottom collisions roughly equals the frequency of 
the bottom-plate, i.e. T < 4.5. Varying the cell width from 100 to 200 bead diameters, no 
influence was reported either. Thus an empirical dispersion relation was proposed to fit all 
the data ranging from H = 3 to H = 9 layers, i.e. A = y/H(X*(d) + g*/f 2 ), with X*(d) = 
7.2 mm, and g* = 1.05 m/s 2 . In contrast, the simulations show that the wavelength does 
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depend on T. In Fig. ^]a, we plot the wavelength L x of a system with H = 6 at constant 
frequency / = 10Hz and for accelerations in the range 2.6 < T < 4.3. 
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FIG. 2. (a) L x as a function of T for H = 6, f = 10Hz, and L/d = 100 (squares), or L/d = 
166 (solid triangles). The error-bars denote the standard deviation of selected averages and the 
horizontal line is the empirical fit of Ref. |]l^] in Figs. 0(a), (b), and (c). the error-bars denote 
the standard deviation of selected averages, (b) L x as a function of / for H = 6, and V = 3.6, L/d 
= 100 (squares), T = 3.6, L/d = 200 (triangles), T = 3.2, L/d = 200 (stars), and T = 3.6, L/d = 
200, H = 10 (x), and H = 14 (circles) . The solid curve is the empirical fit of Ref. |l2| for H = 6. 
(c) L x as a function of H for / = 10Hz, T = 3.6, and L/d = 100 (diamonds) or L/d = 200 (solid 
diamonds), (d) Collection of all simulations from Figs. |2|(b) and (c), in dimensionless scale. The 
dotted line is the dispersion relation expected for gravity waves. 



The solid flat line is the empirical estimate of Ref. [I2|. For L/d = 100 (squares) a 
modulation of the wavelength is clearly evidenced and also for L increased by a non-integer 
factor, i.e. L/d = 166 (solid triangles), this modulation exists, however with a different 
structure. We give some typical standard deviations of the averages, to show that the 
structure is not only due to noise. Since L x is averaged over many periods, not much 
quantitative can be said about L x {t). For the simulations with a large standard deviation 
we mostly observe a significant periodic change of L x {t) from period to period between two 
values. In Fig. |2|b we plot L x versus / for various simulations with H = 6, 10, 14, T = 3.2, 3.6 
and L/d 



100, 200. Comparing the results with the empirical fit of Ref. |12[ for H = 6, 
we observe qualitative agreement with the corresponding simulations. Now, we vary the 
height of the layer from H = 5 to H = 15. The acceleration is kept constant at T — 3.6 
and the system width is L/d = 100 (diamonds) or L/d = 200 (solid diamonds). In Fig. 
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|2]c, the wavelength is plotted as a function of H . We observe an increase of L x with H but 
we cannot extract a functional behavior from our data due to the strong fluctuations. In 
Fig. 0d, the dispersion relation for fluid gravity waves is tested fll3"f . We plot uj 2 /(4gk) as 



a function of hk with h = \/3dH/2, and the wavenumber k = 2n/L x , for the simulations 
of Figs. 0b and c. The dotted line is the expected dispersion relation for gravity waves 
co 2 /(4gk) = ta.nh(hk). Even when the data seem to gather near the line, the fluctuations 
are too strong to allow for a conclusive statement on the dispersion relation. 

In this letter, we present simulations of vibrated 2D layers of grains, using an algorithm 
based on an event driven procedure with a restitution coefficient depending on the impact 
velocity and a collision frequency dependent dissipation threshold. The simulation proce- 
dure keeps us, by algorithmic construction, out of the regime where multiparticle effects 
become dominant. However, we verify that the patterns also occur with a standard inter- 
action model. The present procedure is designed to explore more efficiently larger domains 
of parameter space, i.e. the divergence of the collision frequency (the "inelastic collapse") 
can be avoided. The reported patterns are in qualitative agreement with the experimental 
findings. We observe standing peak patterns at the layer top and - depending on different 
parameters - an arch structure forming at the bottom. The patterns oscillate with twice 
the bottom plate period. The pattern wavelength was systematically extracted using the 
horizontal density correlation function. Like in the experiments, we evidence a regime where 
the wavelength decreases when the frequency increases and an almost constant wavelength 
in the limit of large frequencies. We observe a modulation of the wavelength around the ex- 
perimental empirical determinations, which is triggered by a resonant effect between the box 
size and the bottom plate acceleration. This effect was not reported in previous experiments 



in large cells but something analogous was reported for small cells |13j. For layer heights 
between H = 5 and H = 15, we measure a weak increase of the wavelength compatible with 
previous experimental determinations but the data does not allow for a definite quantitative 
conclusion. 

Finally, we remark that the picture of a completely inelastic block - often used to describe 
a dissipative thick granular layer - is not valid if the system succeeds to choose a state, 
i.e. standing waves, in which energy is not totally dissipated during the contact with the 
bottom. The aim of fully understanding the instability presented here is to unravel the 
physics governing the modes of transport of mass, momentum and energy in a vibrated 
granular material. An open and challenging question is to extract from these parametric 
excitation studies what might be specific to granular assemblies and what can be understood 
in the general framework of hydrodynamic instabilities. The qualitative convergence we find 
here between experimental results and numerical computations is encouraging and calls for 
more detailed studies on both sides. 
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